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ABSTRACT 

We compute the expected low-frequency gravitational wave signal from coalescing massive black hole 
(MBH) binaries at the center of galaxies in a hierarchical structure formation scenario in which seed 
holes of intermediate mass form far up in the dark halo "merger tree". The merger history of dark 
matter halos and associated MBHs is followed via cosmological Monte Carlo realizations of the merger 
hierarchy from redshift z = 20 to the present in a ACDM cosmology. MBHs get incorporated through 
halo mergers into larger and larger structures, sink to the center owing to dynamical friction against the 
dark matter background, accrete cold material in the merger remnant, and form MBH binary systems. 
Stellar dynamical (three-body) interactions cause the hardening of the binary at large separations, while 
gravitational wave emission takes over at small radii and leads to the final coalescence of the pair. A 
simple scheme is applied in which the "loss cone" is constantly refilled and a constant stellar density core 
forms because of the ejection of stars by the shrinking binary. The integrated emission from inspiraling 
MBH binaries at all redshifts is computed in the quadrupole approximation, and results in a gravitational 
wave background (GWB) with a well defined shape that reflects the different mechanisms driving the 
late orbital evolution. The characteristic strain spectrum has the standard h c (f) oc f~ 2 ^ 3 behavior only 
in the range / = 1CV 9 — 10~ 6 Hz. At lower frequencies the orbital decay of MBH binaries is driven 
by the ejection of background stars ("gravitational slingshot"), and the strain amplitude increases with 
frequency, h c (f) oc /. In this range the GWB is dominated by 10 9 - 10 10 M o MBH pairs coalescing 
at < z < 2. At higher frequencies, / > 10~ 6 Hz, the strain amplitude, as steep as h c (f) oc / _1 ' 3 , 
is shaped by the convolution of last stable circular orbit emission by lighter binaries (10 2 — 1O 7 M ) 
populating galaxy halos at all redshifts. We discuss the observability of inspiraling MBH binaries by a 
low-frequency gravitational wave experiment such as the planned Laser Interferometer Space Antenna 
(LISA). Over a 3-year observing period LISA should resolve this GWB into discrete sources, detecting 

60 (w 250) individual events above a S/N = 5 (S/N = 1) confidence level. 

Subject headings: black hole physics - cosmology: theory - early universe - general relativity - 
gravitational waves 



1. INTRODUCTION 

Studies of gravitational wave (GW) emission and of its 
detectability are becoming increasingly topical in astro- 
physics. Technological developements of bars and inter- 
ferometers bring the promise of a future direct observa- 
tion of gravitational radiation, allowing to test one of the 
most fascinating prediction of General Relativity, and, at 
the same time, providing a new powerful tool in the as- 
tronomical investigation of highly relativistic catastrophic 
events, such as the merging of compact binary systems and 
the collapse of massive stellar cores. 

Massive black hole (MBH) binaries are among the pri- 
mary candidate sources of GWs at mHz frequencies (see, 
e.g., Haehnelt 1994; Jaffe & Backer 2003; Wyithe & Loeb 
2003), the range probed by the space-based Laser Inter- 
ferometer Space Antenna (LISA). Today, MBHs are ubiq- 
uitous in the nuclei of nearby galaxies (see, e.g., Magorrian 
et al. 1998). If MBHs were also common in the past (as 
implied by the notion that many distant galaxies harbor 
active nuclei for a short period of their life), and if their 
host galaxies experience multiple mergers during their life- 
time, as dictated by popular cold dark matter (CDM) hi- 
erarchical cosmologies, then MBH binaries will inevitably 



form in large numbers during cosmic history. MBH pairs 
that are able to coalesce in less than a Hubble time will 
give origin to the loudest GW events in the universe. 

Provided wide MBH binaries do not "stall" , their GW 
driven inspiral will follow the merger of galaxies and pre- 
galactic structures at high redshifts. A low-frequency de- 
tector like LISA will be sensitive to waves from coalescing 
binaries with total masses in the range 10 3 — 10 6 Mq out 
to z ~ 5 — 10 (Hughes 2002). Two obvious outstanding 
questions are then how far up in the dark halo merger 
hierarchy do MBHs form, and whether stellar dynamical 
processes can efficiently drive wide MBH binaries to the 
GW emission stage. The luminous z « 6 quasars dis- 
covered in the Sloan Digital Sky Survey (Fan et al. 2001) 
imply that black holes more massive than a few billion so- 
lar masses were already assembled when the universe was 
less than a billion years old. These MBHs could arise nat- 
urally from the growth by Eddington-limited gas accretion 
of seed stellar-mass or intermediate-mass holes forming at 
z > 10 (Haiman &; Loeb 2001). It seems also clear that, 
following the merger of two halo+MBH systems of com- 
parable mass ("major mergers"), dynamical friction will 
drag in the "satellite" halo (and its MBH) toward the cen- 
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ter of the more massive progenitor: this will lead to the 
formation of a bound MBH binary in the violently relaxed 
core of the newly merged stellar system. What remains 
uncertain is the demography of the MBH population as a 
function of redshift and the dynamics of MBH binaries, i.e. 
whether binaries will "hung up" before the backreaction 
from GW emission becomes important. 

In this paper we study the expected gravitational wave 
signal from inspiraling binaries in a hierarchical structure 
formation scenario in which seed holes of intermediate 
mass form far up in the dark halo "merger tree" . The 
model has been discussed in details in Volonteri, Haardt 
& Madau (2003, hereafter Paper I). Seed holes are placed 
within rare high-density regions (minihalos) above the cos- 
mological Jeans and cooling mass at redshift 20. Their 
evolution and growth is followed through Monte Carlo 
realizations of the halo merger hierarchy combined with 
semi-analytical descriptions of the main dynamical pro- 
cesses, such as dynamical friction against the dark mat- 
ter background, the shrinking of MBH binaries via three- 
body interactions, their coalescence due to the emission of 
gravitational waves, and the "gravitational rocket" . Ma- 
jor mergers lead to MBH fueling. The long dynamical fric- 
tional timescales leave many MBHs "wandering" in galaxy 
halos after a minor merger. Bound pairs form after major 
mergers and lose orbital angular momentum by capturing 
stars passing within a distance of the order of the binary 
semi-major axis and ejecting them at much higher veloc- 
ities. The heating of the surrounding stars by a decaying 
MBH pair creates a low-density core out of a preexisting 
stellar cusp, slowing down further binary hardening (see, 
e.g., Milosavljevic & Merritt 2001). The model reproduces 
rather well the observed luminosity function of optically- 
selected quasars in the redshift range 1 < z < 5 (Paper 
I), and provides a quantitative explanation to the stel- 
lar cores observed today in bright ellipticals as a result of 
the cumulative eroding action of shrinking MBH binaries 
(Volonteri, Madau, & Haardt 2003, hereafter Paper II). 

The paper is organized as follows. In § 2 we describe 
our merger tree algorithm coupled with semi-analytical 
recipes designed to track the dynamical history of MBHs 
in a cosmological framework. We compute black hole bi- 
nary coalescence rates as a function of redshift and mass, 
and show that MBH binaries will shrink to the gravita- 
tional wave emission stage by scattering off background 
stars from stellar cusps. In § 3 we summarize the basic 
theory of GW emission. In § 4 we compute the cosmolog- 
ical GW background (GWB) from inspiraling MBH bina- 
ries and discuss the observability of individual coalescence 
events by LISA. Finally, we summarize our main results 
in § 5. Unless otherwise stated, all results shown below 
refer to the currently favoured ACDM world model with 
VL M = 0.3, n A = 0.7, h = 0.7, n b = 0.045, a 8 = 0.93, and 
n = 1. 

2. MBH BINARIES 

2.1. Assembly and growth of MBHs 

Following Papers I and II, we track backwards the 
merger history of parent halos with present-day masses 
in the range 10 11 - 10 15 M with a Monte Carlo al- 
gorithm based on the extended Press-Schechter formal- 
ism (see, e.g., Cole et al. 2000). The tree is averaged 



over a large number of different realizations, for a total 
of 220 simulated halos, and the results weighted by the 
Press-Schechter halo mass function at the chosen output 
redshift. Prcgalactic seed holes form with intermediate 
masses (to, = 150 M Q ) as remnants of the first genera- 
tion of massive metal-free stars with to* > 260 M Q that 
do not disappear as pair-instability supernovae (Madau & 
Rees 2001). We place them in isolation within halos above 
Mgeed = 1-6 x 10 7 M Q collapsing at z = 20 (the highest 
redshift computational costs allow us to follow the merger 
hierarchy to) from rare > 3.5<r peaks of the primordial 
density field. While Z = stars with 40 < to* < 140 M 
are also predicted to collapse to MBHs with masses ex- 
ceeding half of the initial stellar mass (Heger & Woosley 
2002), we find that the GWB signal in the LISA window 
is not very sensitive to the precise choice for the seed hole 
mass. 

Hydrodynamic simulations of major mergers have shown 
that a significant fraction of the gas in interacting galaxies 
falls to the center of the merged system (Mihos & Hern- 
quist 1996): the cold gas may be eventually driven into 
the very inner regions, fueling an accretion episode and 
the growth of the nuclear MBH (see, e.g. Kauffmann & 
Haehnelt 2000). Since the local MBH mass density is con- 
sistent with the integrated luminosity density of quasars 
(Yu & Tremaine 2002), the fraction of cold gas ending up 
in the hole must depend on the properties of the host halo 
in such a way to ultimately lead to the observed correlation 
between stellar velocity dispersion and SMBH mass. As in 
Papers I and II, we assume that in each major merger (de- 
fined here as a merger between halos of mass ratio > 0.3), 
the more massive hole accretes, at the Eddington rate, a 
gas mass that scales with the fifth power of the circular 
velocity of the host halo, 

ATO acc = 3.6 x 10 6 M AC V c % a , (1) 

where V c ,i50 is the circular velocity of the merged system 
in units of 150 kins" 1 . This relation follows from com- 
bining the tobh — and the V c — cr* relations given in 
Fcrrarcsc (2002, see Papers I and II). The normalization 
factor JC is of order unity and is fixed a posteriori in or- 
der to reproduce the correlation between stellar velocity 
dispersion and nuclear black hole mass observed today in 
nearby galaxies (Ferrarese & Merritt 2000; Gebhardt et al. 
2000). 

Inside a halo of mass M, the dark matter is distributed 
according to a NFW profile (Navarro, Frenk, & White 
1997), 

. , M , n . 

PDM(r) = f— -, 2 . r , (2) 

47rr(r + r vir / c) z J (c) 

where r v ; r is the virial radius, c is the halo concentration 
parameter and /(c) = ln(l + c) — c/(l + c). The distribu- 
tion of concentrations and its scaling with halo mass and 
redshift of collapse are taken from Bullock et al. (2001). 
Dynamical friction drags satellite progenitors - and their 
MBHs if they host one or more - towards the center of 
the more massive preexisiting system, on a timescale that 
depends on the orbital parameters of the satellites (van 
den Bosch et al. 1999) and on tidal mass loss/evaporation 
(Taffoni et al. 2003). The initial occupation fraction of 
seed holes is so small that the merging of two halos both 
hosting a black hole is a rare event even with the assumed 
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"bias" , and MBHs evolve largely in isolation. At redshifts 
below 15, however, more than 10% of hosts contain two or 
more MBHs, and bound binaries start forming in signifi- 
cant numbers. 

The evolution of a bound binary is determined by the 
initial central stellar distribution. We model this as a 
(truncated) singular isothermal sphere (SIS) with one- 
dimensional velocity dispersion cr* and density 



P*{r) 



2irGr 2 ' 



(3) 



The stellar velocity dispersion is related to the halo cir- 
cular velocity V c at the virial radius following Ferrarese 
(2002), 

logV c = (0.88 ±0.17) log cr* + (0.47 ±0.35). (4) 

2.2. Binary merger rates 

The scmimajor axis a(t) of a bound binary continues 
to shrink owing to dynamical friction from distant stars 
acting on each MBH individually, until it becomes "hard" 
when the separation falls below (Quinlan 1996) 



Gm 2 _ 2 
a h = = 0.2 pc m 6 <r 70 , 



(5) 



where ni2(< mi) is the mass of the lighter hole, mg = 
m 2 /10 6 M©, and (T70 is the stellar velocity dispersion in 
units of 70 kms -1 . After this stage the MBH pair hardens 
via three-body interactions, i.e., by capturing and ejecting 
at much higher velocities the stars passing by within a 
distance <~ a (Bcgelman, Blandford, & Rees 1980). We 
assume that the "bottleneck" stages of binary shrinking 
occur for separations a < ah,; during a galactic merger, 
after a dynamical friction timescale, we place the MBH 
pair at and let it evolve. The hardening of the binary 
modifies the stellar density profile, removing mass interior 
to the binary orbit, depleting the galaxy core of stars, and 
slowing down further hardening. If M c j is the stellar mass 
ejected by the pair, the binary evolution and its effect on 
the galaxy core are determined by two dimcnsionless quan- 
tities: the hardening rate 



H = 



cr* d 1 
Gp* dt a ' 



and the mass ejection rate 



J 



1 



dM c 



(mi + m 2 ) d\n(l/a) ' 



(6) 



(7) 



The quantities H and J can be found from scattering ex- 
periments that treat the test particle-binary encounters 
one at a time (Quinlan 1996). We assume as in Papers I 
and II (see also Merritt 2000) that the stellar mass re- 
moval creates a core of radius r c and constant density 
pc = P*(r c ), so that the total mass ejected as the binary 
shrinks from ah to a can be written as 



M cj = -^-(r c - n) + Mi ~M C = - 



4^(r 



G 



(8) 



where M c = 4irp c r*/3, M t = ATrpirf/3, n = r c {t = 0) 
is the radius of the (preexisting) core when the hardening 



phase starts at t = 0, and pi = p*(r^). The core radius 
then grows as 



3 „, , f ah J{a) 



r c (t) = r t + —G(m x + m 2 ) 



f ah J(a 

Ja(t) a 



da. (9) 



The binary separation quickly falls below r c and subse- 
quent evolution is slowed down due to the declining stellar 
density, with a hardening time, 



th = \a/a\ 



2^r c (tf 
Ha* a ' 



(10) 



that becomes increasingly long as the binary shrinks. The 
mass ejected increases approximately logarithmically with 
time, and the binary "heats" background stars at radii 
r c 3> a. We assume that the effect of the hierarchy of 
MBH binary interactions is cumulative (the "core preser- 
vation" model of Paper II), i.e. the mass displaced by the 
pair is not replenished after every major merger event and 
j-j continues to grow during the cosmic evolution of the 
host. This is supported by N-body simulations involving 
mergers of spherical galaxy models with different density 
profiles, and showing that the remnant profile is quite close 
to the profile of the progenitors - in other words that the 
core appears to be preserved during such mergers (e.g. Ful- 
ton & Barnes 2001). The above scheme can be modified 
by rare triple black hole interactions, when another major 
merger takes place before the pre-existing binary has had 
time to coalesce (see, e.g., Hut & Rees 1992). In this case 
there is a net energy exchange between the binary and the 
third incoming black hole, resulting in the ejection of the 
lighter hole and the recoil of the binary. The binary also 
becomes more tightly bound. The reader is referred to Pa- 
per I for a detailed discussion of the role played by triple 
interactions. If hardening continues down to a separation 



a gw = 0.0014 pc 



(mi + m 2 )mim 2 



10 18 - 3 M Q 3 



1/4 



f V4 



(11) 



the binary will coalesce within tg Gyr due to the emis- 
sion of gravitational waves. An equal mass pair must then 
manage to shrink by a factor 



ah 



J-gW 



150 m^ 4 er 70 2 t g 



1/4 



(12) 



for gravity wave emission to become efficient. Since the 
hardening and mass ejection rate coefficients arc H w 15 
and J w 1 in the limit of a very hard binary (Quinlan 
1996), the hardening time in equation (10) can be rewrit- 
ten (for mi = m 2 and n — ► 0) as 



t h ~ 4.7 x 10 4 yr m 6 cr 70 3 {a h /a) \n{a h /a) 



5.6 x 10 4 yr m°' 3D (a h /a) ln z (a h /a), 



(13) 



where the second equality assumes the wbh — relation of 
Ferrarese (2002). Allowing for the cumulative effect of a hi- 
erarchy of MBH binary interactions (i.e. for the existence 
of a pre-existing core of radius rj) increases somewhat the 
hardening time in equation (13). Yet, we find th to be com- 
parable or shorter than the then Hubble time at all but the 
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highest redshifts. Figure 1 shows the number of MBH bi- 
nary coalescences per unit redshift per unit observed year, 
dN/dzdt, predicted by our model. The observed event rate 
is obtained by dividing the rate per unit proper time by 
the (1 + z) cosmological time dilation factor. Each panel 
shows the rate for different tobh = "*i + Tt2 mass inter- 
vals, and lists the integrated event rate, dN/dt, across the 
entire sky. The number of events per observed year per 
unit redshift peaks at z = 2 for 10 7 < tobh < 10 9 M Q , 
at z = 3 - 4 for 10 5 < m BH < 10 7 M Q , and at z = 10 
for ?tjbh < 10 5 M©, i.e. the lower the black hole mass, 
the higher the peak redshift. Beyond the peak, the event 
rate decreases steeply with cosmic time. Previous calcu- 
lations of this quantity have either assumed instantaneous 
black hole coalescence and computed halo merger rates 
only out to z <~ 5 (Menou, Haiman, & Narayanan 2001), 
or have assigned a constant (independent of black hole 
mass, background stellar density, and redshift) efficiency 
factor to the coalescence process (Wyithe & Loeb 2003), 
or have used phenomenological and empirical rates (Jaffe 
& Backer 2003). Our rates are significantly smaller than 
those computed by Wyithe & Loeb (2003), consistent with 
their assumption that all MBH binaries actually coalesce. 
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sity that is allowed to decrease following equation (9), not 
the density of low angular momentum stars. The effect of 
loss-cone depletion (the depletion of low-angular momen- 
tum stars that get close enough to extract energy from a 
hard binary) has traditionally been one of the major uncer- 
tainties in constructing viable merger scenarios for MBH 
binaries (Begelman et al. 1980). Several processes have 
been recently studied that should all mitigate the prob- 
lems associated with loss-cone depletion, such as the wan- 
dering of the binary center of mass from the galaxy center 
induced by continuous interactions with background stars 
(Chatterjee, Hernquist, & Loeb 2003), the large supply of 
low-angular momentum stars in significantly flattened or 
triaxial galaxies (Yu 2002), the presence of a third hole 
(Blaes, Lee, & Socrates 2002), and the randomization of 
stellar orbits due to the infall of small satellites (Zhao, 
Hachnelt, & Rees 2001). Our scheme also neglects the 
role of gaseous - rather than stellar dynamical - processes 
in driving the evolution of a MBH binary (see, e.g., Escala 
et al. 2004). We will try to address the effect of some these 
uncertainties on the predicted GWB in § 4. 

There is another process that our hierarchical MBH evo- 
lution scenario includes: this is the "gravitational rocket" , 
the recoil due to the non-zero net linear momentum car- 
ried away by GWs in the coalescence of two unequal mass 
black holes. Radiation recoil is a strong field effect that 
depends on the lack of symmetry in the system, and may 
diplacc MBHs from galaxy centers or even eject them into 
intergalactic space (Redmount & Rees 1989). Its outcome 
is still uncertain, as fully general relativistic numerical 
computations of radiation reaction effects are not avail- 
able at the moment. Here we have used the extrapolated 
perturbative results on test mass particles of Fitchett & 
Detwcilcr (1984), predicting in a Schwarzschild geometry 
recoil velocities larger than 200 km s" 1 only for mass ratios 
m2/mi > 0.4. As shown in Figure 2, in our model the typ- 
ical mass ratio of merging binaries exceeds 0.4 at z > 10. It 
is at early times then that the gravitational rocket is effec- 
tive at ejecting MBHs from the shallow potential wells (es- 
cape velocities less than 100 kms -1 ) of their hosts (Madau 
et al. 2004). Overall, we find that w 25% of coalesced pairs 
recoil and escape their host halos at z ~ 10, and that the 
ejected holes have typical masses of > 1000 M Q . The de- 
tectability of recoiling MBHs has been recently addressed 
by Madau & Quataert (2004). 



3. GRAVITATIONAL WAVE SIGNAL 



Fig. 1. — Number of MBH binary coalescences observed per year 
at z = 0, per unit redshift, in different tobh = mi + TO2 mass 
intervals. Each panel also lists the integrated event rate, dN/dt, 
predicted by our model. The rates (solid lines) are compared to a 
case in which triple black hole interactions are switched off (dotted 
lines). Triple hole interactions increase the coalescence rate at very 
high redshifts, while, for 10 < z < 15, the rate is decreased because 
of the reduced number of surviving binaries. Dashed lines: rates 
computed assuming binary hardening is instantaneous, i.e. MBHs 
coalesce after a dynamical friction timescale. 

While our calculations are the first to include a de- 
tailed model of MBH binaries dynamics, we note that our 
scheme implicitly neglects the depopulation of the "loss 
cone" (Frank & Rees 1976), since it is the total stellar den- 



3.1. Basic theory 

The theory of GW emission in the quadrupole approxi- 
mation is textbook topic (e.g., Schutz & Ricci 2001; Mag- 
giore 2000), and we recall here only the basic relations. 
The Einstein equations can be linearized over the per- 
turbed Minkowski metric. The perturbed part of the met- 
ric or "strain" admits a wave solution out of the source of 
the stress-energy tensor. Such solution becomes particu- 
lary easy to treat in the so-called transverse-traceless (TT) 
gauge. In the TT gauge the strain at distance r from the 
source has the form 



h T J = 



(t,r) 



2G 



Hk 



(* 



(14) 



5 



so that only the spatial components are not vanishing. The 
term QjiJ is the reduced quadrupole momentum projected 
into the TT gauge, containing the stress-energy tensor of 
the gravitational source, and evaluated at the retarded 
time t — r/c. 
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Fig. 2. — Normalized distribution of mass ratios of coalescing 
MBH binaries at different epochs, averaged over 120 realizations. 
At very high redshift binary members are seed holes with nearly 
equal masses; as time goes on MBHs grow due to gas accretion, and 
low mass ratios become more probable. 

It is then possible to associate an energy momentum 
pseudo-tensor to the wave, which, to first perturbative or- 
der, transforms as a tensor. Such pseudo-tensor is 



t 



32ttG 



(15) 



where <> is an average over wavelength. If we consider a 
GW propagating along the positive x direction, then the 
associated energy flux will be 



dE, 



TT h TTij\ 



(16) 



Taking equation (14), deprojecting the quadrupole mo- 
mentum Q TT out of the TT gauge to get Q, and integrat- 
ing over the surface S, one derives the luminosity associ- 
ated with GW emission: 



dE s 



dt 



dE, 



gw 



dtdS 



dS 



G_ 

5c 5 



d 3 



dt 3 ^ J dt 3 



d 3 



(17) 



This depends on the time third derivative of the (reduced) 
quadrupole momentum of the source. 

3.2. Gravitational waves from binary systems 

In the stationary case, i.e., assuming no orbital decay, 
the GW emission spectrum of a MBH pair in a circular or- 
bit of radius a is a delta function at rest-frame frequency 



f r = ui/iz, where ui = \J G{m\ + m^J/ 'a 3 is the Keplerian 
angular frequency of the binary. Orbital decay (due to GW 
emission and/or the ejection of background stars) results 
in a shift of the emitted frequency to increasingly larger 
values as the binary evolution proceeds. The energy spec- 
trum integrated over the entire radiating lifetime of the 
source (a continuum formed by the superposition of delta 
functions) can be computed once the energy emitted per 
logarithmic frequency interval, 



dE z 



dE gw dt da 

'fr 



(18) 



rfln/ r dt dadf r r ' 

is known. This is related to the characteristic strain, h. 
through 



dE, 



gw 



d\nf r 
where h r is defined as 



*f?h 2 c (f r ), 



(h af3 (t)h a "(t)) = 2 (dlnf r )h 2 c (f r ). 



(19) 



(20) 



The characteristic strain is related to the Fourier trans- 
form of the strain h(f r ), h 2 c {f r ) = 2f r \h(f r )\ 2 . In the case 
of a binary system equation (17) can be rewritten as: 



dE, 



gw 



dt 



32^10/3 G 7/3 

5c 5 



(Mf r ) 



10/3 



(21) 



The radiated power is an increasing function of frequency 
and of the "chirp mass" of the system, 



M = 



in 



3/5 3/5 



(mi + m 2 ) 



1/5 1 



(22) 



The term da/df in equation (18) is readily obtained by 
differentiating the relation f r — lo/tt, 



da 
df r 



G(mi 



m 2 J 



1/3 



-5/3 



(23) 



Finally, the term da/dt represents the orbital decay due to 
angular momentum losses. When the backreaction from 
GW emission dominates, the orbit shrinks at a rate 



da\ 
dt) i 



64 G 3 



5 c 5 a 3 

and the resulting spectrum of GWs is 



m\mi{m\ + ni2), 



dE gw 
d In f r 



(ttG) 



2/3 



M 5/3 f 2/3 



(24) 



(25) 



gw 



In the early stages of binary evolution, however, angu- 
lar momentum losses are driven by stellar dynamical pro- 
cesses, and the emerging GW spectrum will show a dif- 
ferent scaling with frequency. The model outlined in the 
previous section leads to the following expression for the 
orbital decay rate in the "gravitational slingshot" regime 
(eq. 10): 

Ha*a 2 (2g) 



dt 



27T7-2 
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Assuming an initially cuspy profile (r^ = 0) and in the 
limit of a very hard binary (J w 1), the core radius (eq. 
9) created by the shrinking pair can be written as 



r c (t) « -^G(mi +m 2 ) \n(a h /a) 



(27) ^ 



and the hardening rate becomes proportional (neglecting 
the logarithmic dependence) to 



dd \ c , x 9 9 

— \ cx <(mi + m 2 ) a 



(28) 



Contrary to the case when GW losses dominates, the bi- 
nary decay slows down with time because of the declining 
stellar density, and the rate \da/dt\ decreases as the orbit 
shrinks. Inserting the above expression into equation (18) 
yields the following scaling for the emitted GW spectrum 
in the stellar ejection regime: 



dEg W 



(x m 1 m 2 {m 1 + m 2 ) a„ f r . 



(29) 



The GW spectrum arising from a binary system during 
the slow inspiral phase spans a broad but finite interval of 
emitted frequencies, / m j n < f r < /max- We assume that 
the binary starts emitting GWs at separation a = ah, so 
that the lower limit / m i n is 



Q 3 

/mi„ = ^(mi + m 2 ) 1/2 m 2 3/2 . (30) 



The upper limit / ma x is set by the frequency emitted at 
a = 6Gmi/c 2 , the radius of the closest stable circular orbit 
for two non-rotating holes, 



/max — 



6 3 / 2 ttG 



(mi + m 2 ) 1/2 m 1 3/<2 . 



(31) 



We do not address in this work the higher frequency radi- 
ation emitted during the subsequent coalescence and ring- 
down phases (Flanagan & Hughes 1998). Figure 3 shows 
three examples of GW spectra emitted by inspiraling bina- 
ries. The two regimes, stellar slingshot at low frequencies 
and GW losses at high frequencies, are clearly recogniz- 
able. 




-10 -5 
log emitted frequency [Hz] 

Fig. 3. — Energy per natural log frequency interval around 
f r emitted in GWs by an equal mass (mi = IH2) black 
hole binary system in a stellar background of velocity disper- 
sion (j * . From top to bottom, the three curves corresponds 
to (mi, 0-,,) = (10 8 M Q , 180 kms" 1 ), (10 6 M ,66 kms" 1 ), and 
(10 4 M Q ,24 kms -1 ), respectively. 



3.3. Total background from cosmic sources 

We can now compute the total background - random 
GWs arising from a large number of independent events - 
produced by a population of sources at cosmological dis- 
tances. Following Phinncy (2001), the present-day energy 
density in GWs of frequency / per natural log frequency 
interval is 



c 2 p gw (/) - 



! hl(f) = f 
Jo 



dz 



N(z) dE s 



1 + z d\nf r 



(32) 



where N(z)dz is the comoving number density of catas- 
trophic events in the redshift interval z, z + dz, the factor 
1/(1 + z) accounts for the redshifting of gravitons, and 
the observed frequency / is related to the rest-frame fre- 
quency f r by / = / r /(l + z). Equation (32) assumes that 
all sources emit the same intrinsic spectrum (i.e. all bi- 
naries have the same mass and shrink in the same stellar 
background), and that the inspiraling timescale is much 
shorter than the then Hubble time, i.e. the inspiral oc- 
curs at the same redshift as the coalescence. When binary 
decay occurs on timescales that are comparable or longer 
than the expansion time, equation (32) can be generalized 
as 



c 2 Pgw(/) = / 



(33) 



where N(f r , t) is the comoving number density of binaries 
at cosmic time t that radiate at frequencies between /,. 
and f r + df r , dt/dz = —(1 + z)^ 1 H(z)^ 1 where H(z) is 
the Hubble parameter, and dE gV! /dt is the GW luminosity 
at fr = /(I + z) measured in the rest-frame of the pair. 
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The density of sources N(f r ) is related to the space den- 
sity of binaries with separations between a and a + da by 
N(f r )df r = -N(a)da. If N b is the birth rate (after dy- 
namical friction) of MBH binaries, and N c is the rate of 
coalescences, then N(a) is given by the continuity equation 



dN 
~dt 



djaN) 
da 



N b 5(a-a h )-N c 5(a). (34) 



Note that, since in the situation we are considering the 
sources of GWs are MBH binaries with a mass function 
that changes with cosmic epoch, the term N(f r , t) dE gV! /dt 
must be summed up over all binary types. 



3.4. Resolvable signal 

While the data stream of LISA will include confusion 
noise arising from a large number of unresolved galactic 
and extragalactic white dwarf binaries (Farmer & Phinney 
2003; Nelemans, Yungelson, & Portegies-Zwart 2001), the 
GWB from MBH binaries will be resolvable into discrete 
sources at least in the high frequency end of the LISA sen- 
sitivity range. The nature (stochastic or resolved) of the 
GWB can be assessed by counting the number of events 
per resolution frequency bin whose observed signal is larger 
than a given sensitivity threshold /i c ,min- The background 
is resolved if such number is less than unity and the cumu- 
lative signal from all sources below threshold is negligible. 

The resolution frequency bin is set by the minimum fre- 
quency resolvable by a given mission, and it is simply the 
inverse of the mission lifetime. The characteristic strain 
h c of a source at comoving coordinate distance r(z) is 



h\/n, 



(35) 



where the strain amplitude h (sky-and-polarization aver- 
aged) is given by 



8tt 2 / 3 G^M 5 '* 2/3 
~ 10V2 c 4 r(z) tr ' 



(36) 




-6 -4 -2 
log observed frequency [Hz] 

Fig. 4. — Gravitational wave background from inspiraling MBH 
binaries in the hierarchical scenario described in the text. The 
square of the characteristic strain (directly observable by a long- 
base interferometer) il plotted vs. wave frequency (solid line). Thick 
solid line: LISA single-arm Michelson sensitivity curve. Thick dash 
at / ~ 10~ 9 : current limits from pulsar timing experiments. Long- 
dashed line at 10~ 6 < / < 0.1 Hz: expected strain from extra- 
galactic white dwarf binaries (Farmer & Phinney 2003). Dotted line 
at 10~ 4 < / < 10~ 3 Hz: expected strain from unresolved galactic 
white dwarf binaries (Nelemans et al. 2001). 



and n is the number of cycles a source spends at frequency 
f r , i.e., n = f 2 /f r (Thorne 1987). 3 Note that, for a fi- 
nite observation time r, the number of cycles actually ob- 
served at frequency f r can not exceed fr. The frequency 
evolution f r can be obtained by combining equations (23) 
and (24) in the GW regime, (23) and (26) in the slingshot 
regime. At the endpoint of the inspiraling phase, the num- 
ber of cycles spent near that frequency is close to unity, so 
that h c ~ h. 

3 The formula for h in Thorne (1987) contains a multiplicative extra factor 4/3 to account for the Earth rotation. 



4. RESULTS 

4.1. Gravitational Wave Background 

Merger tree realizations allow us to follow the time evo- 
lution of individual inspiraling binaries, and solve equa- 
tion (33) by performing a direct weighted sum of the 
emission from individual evolving pairs. The resulting 
GWB is shown in Figure 4, where the square of the 
characteristic strain is plotted against observed frequency. 
The peak occurs for / poa k ~ 10 -10 Hz. For / < / poa k 
the spectrum is shaped by the evolution of binaries in 
the stellar slingshot regime, and, according to equation 
(29), hi oc f 2 just below / poa k- At even lower frequen- 
cies the energy distribution is somewhat steeper than f 2 
because of the assumed cutoff in the individual emis- 
sion spectra for a = ah (see Figure 3). For frequen- 
cies / > /peak, but below 10~ 6 Hz, h 2 oc /~ 4 / 3 , as ex- 
pected in the classical GW regime (see eqs. 19 and 25). 
At even higher frequencies, / > 10~ 6 Hz, the strain is 
shaped by the convolution of the sharp frequency cutoff 
at the last stable orbit of individual binaries (see Figure 
3), and in this range one has approximatively h 2 oc J -2 ' 6 . 
The slope in this regime is directly related to the rate 
of binary coalescences and to the MBH mass function. 
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Fig. 6. — GWB from different redshift intervals. From top to 
bottom, the spectrum of the characteristic strain is computed inte- 
grating the emission from events at < z < 1 (solid line), 1 < z < 2 
(long- dashed line), 2 < z < 3 (short- dashed line), and z > 3 (dotted 
line), respectively. 

Lighter binaries show up at higher frequencies, where 
the background drops in amplitude. For / > 10 -4 Hz, 
the signal is dominated by mi < 10 7 M Q binaries, with 
comparable contributions from events at z < 3 and z > 3. 
The peak frequency / pea k of each curve in Figure 6 is de- 
termined by two factors acting in opposite directions, i.e. 
the Hubble expansion and the time evolution of the binary 
mass function. The latter remains approximately constant 
in the redshift interval < z < 2, so that the spectrum 
originating from sources at 1 < z < 2 appears, because of 
the Hubble expansion, redshifted with respect to the sig- 
nal from < z < 1 events. The peak frequency of events 
at z > 2 becomes bluer because the effect of the reduced 
average binary mass at increasing redshifts dominates over 
the Hubble expansion. 

4.2. Stochastic versus resolved signal 

Merger tree realizations allow us to compute the GW 
signal from individual MBH binaries, and to compare it 
with the LISA instrumental noise. Figure 7 depicts the 
number of systems above the S/N = 5 and S/N = 1 sen- 
sivity thresholds, per resolution frequency bin (A/ = 10~ 8 
Hz assuming a 3-year mission duration), and shows that 
such number is less than unity all over the surveyed fre- 
quency range. For reference, we also plot the total number 
of sources per bin, i.e. the number of sources regardless of 
the strength of their GW signal. Even in this ideal case, 
f° r / ^ 10~ 3 ' 5 Hz, the GWB would be resolved into dis- 
crete events. We have also checked that the signal from 
the rare resolved events accounts for more than the 99% of 

4 The curve is taken from http://www.srl.caltech.edu/~shane/sensitivity, where it is given in terms of the "effective strain" h c g = h c /\/J- 

5 The figures given are obtained neglecting every source of noise other than instrumental. In particular, the stochastic noise from galactic 
WD- WD binaries shown in Figure 4 (Lommen et al. 2003) has not been considered here. 



Fig. 5. — Integrated GWB from inspiraling MBH binaries in dif- 
ferent mass ranges. From top to the bottom (as indicated by the 
labels), the curves show the signal produced by events with 9 < 
log(mi/M ) < 10 8 < log(mi/M©) < 9. ..,2 < log(mi/M©) < 3. 
Here mi is the mass of the heavier black hole. 

Figure 4 also shows the LISA single-arm Michelson sen- 
sitivity curve for a sky-and-polarization averaged signal 
with signal-to-noise S/N = l. 4 It is apparent how LISA 
could detect the GWB above 10~ 5 Hz, while the peak at 
nHz frequencies is close to current limits of pulsar tim- 
ing experiments (Lommen 2002). We have also plotted 
the GWB expected in the LISA window from the cosmo- 
logical population of white dwarf binaries computed by 
Farmer & Phinney (2003), and the GWB from unresolved 
galactic white dwarf binaries (Nclemans et al. 2001). Our 
estimate of the background due to MBH binaries lies well 
above the expected strain from cosmological white dwarf 
binaries. 

The contribution to the GWB from MBH binaries in 
different mass ranges is depicted in Figure 5, while Fig- 
ure 6 shows the signal from different redshift intervals. 
From these figures it is clear how, for / < 10~ 6 Hz, the 
background is dominated by events at moderate redshifts, 
< z < 2, from very massive binaries, mi > 10 8 M Q . 
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the total GWB, i.e. the numerous sources below the sensi- 
tivity threshold provide a negligible background at Earth. 
We conclude that the GWB from cosmological MBH bi- 
naries will be resolved into discrete events all over LISA 
frequency range. 

Integrating over frequencies, the number of MBH binary 
systems observable by LISA with S/N in excess of 5 (1) is 
w 20 (200). 5 These numbers refer to stationary events, 
i.e., binaries at large separation whose GW signal does not 
shift much in frequency, hence longlasting compared to the 
observation time. The number of such events is indepen- 
dent on any (plusible) mission duration. On the contrary, 
the number of observable bursts, i.e., binaries during the 
observation period, clearly depends on the mission life- 
time. In 10 8 sees there will be w 200 bursts (see upper 
left panel in Figure 1), and among them, w 35 (50) are 
detectable, by LISA, with S/N > 5 (S/N > 1). 




-5 -A -3 -2 -1 

log observed frequency (Hz) 

Fig. 7. — Number of MBH binaries per 10~ 8 Hz frequency bin 
contributing to the GWB as received today. Dashed line: all bi- 
naries. Upper solid line: binaries above LISA S/N = 1 sensitivity 
threshold. Lower solid line: binaries above LISA S/N = 5 sensitiv- 
ity threshold. 

The redshift distribution of detectable events is shown 
in Figure 8, in the case of S/N — 5 threshold. The distri- 
bution of bursts is similar to that of coalescences shown in 
Figure 1, lower left panel, indicating that a relevant frac- 
tion of bursts at high redshifts produces a signal above 
threshold. The peak of the redshift distribution of sta- 
tionary events is shifted at lower redshift by two different 
instrumental effects. First, the inclusion of a detection 
threshold selects preferentially binaries at low redshifts. 
Second, as for wide binaries the time spent at a given fre- 
quency can be larger than the observation time, only a 
fraction of the signal (see eqs. 35 and 36) may actually be 
detected in a 10 8 s observation. This fact causes the drop 
of the stationary event counts below z ~ 3 seen in Figure 




redshift 



Fig. 8. — Number of events per unit redshift interval resolved 
by LISA with S/N > 5 in 10 s sees. Thick-solid histogram: total 
number of events in 10 8 sees. Solid histogram: number of stationary 
events. These events are of much longer duration compared to the 
mission lifetime. Dashed histogram: number of bursts in 10 8 sees. 
These events are of short duration compared to the mission lifetime. 



4.3. Test and comparison 

It is fair at this stage to ask how sensitive our predic- 
tions are to the details and uncertainties of the adopted 
scheme for binary formation and evolution. To answer this 
question we have run two test models, arbitrarily dividing 
and multiplyng the computed hardening time th by a fac- 
tor of 3, i.e. significantly increasing or reducing the rate 
of binary decay by stellar dynamical processes. 

In the "fast" hardening case, / pca k increases by a factor 
<~ 30% and the peak amplitude drops by ~ 10%. This 
is due to the fact that, when th is shorter, GW emission 
losses dominates over stellar slingshots only at shorter bi- 
nary separations. The opposite occurs in the case of "slow" 
hardening, where / pea k is lowered by a factor ~ 30 %, and 
the peak amplitude is ~ 20% larger than in the standard 
case. For / > / pea k, the background amplitude is 15% 
larger in the "fast" hardening case, because a larger num- 
ber of MBH binaries can now coalesce. The coalescence 
rate is now 78 per observed year, to be compared to 64. 
Again, the opposite is true in case of "slow" hardening, 
the coalescence rate decreasing to 50 per observed year. 
The situation is different in the LISA window. A faster or 
slower hardening has basically no effects, as the increased 
or reduced rate of coalescence affects very massive bina- 
ries only, which do not contribute to the LISA frequency 
range. Consequently, the precise value of the hardening 
time has no significant impact on the number of sources 
(cither stationary or bursts) observable by LISA. This is, 
of course, only true as long as th is less than the then 
Hubble time. If MBH binaries were instead to "stall" (be- 
cause of the depopulation of the loss cone), the rate of 
detectable events could become negligible. We also find 
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our predictions to be sensitive to the occupation fraction 
of halos hosting nuclear MBHs. Within our assumptions 
(rare "high-cr" seed holes), this is of order 50% at z ~ 5 
for halos more massive than 10 10 M©. A lower occupation 
fraction would decrease the rate of inspiraling events, at 
the expenses of making MBHs non- ubiquitous in the nuclei 
of nearby galaxies. The black hole coalescing rate would 
be higher if the seeds were more common at early times. 

We have compared our results with the calculations of 
Wyithe & Loeb (2003), who estimated the GWB from 
MBH binaries using a different, but similar in spirit, ap- 
proach. In their work the time dependent binary mass 
function is explicitly derived from the extended Press- 
Schechtcr formalism. The nature of such methodology 
hampers the possibility of following the evolution of in- 
dividual MBH pairs (e.g., triple interactions as well as the 
differential frequency broadening due to the non-negligible 
duration of the orbital decay phase are neglected). The 
two estimates of the GWB appear in good agreement (to 
within a factor 2) at high frequencies (/ > 10 -4 Hz), a 
fact reflecting the common hierarchical scenario, the sim- 
ilar normalization to the present-day tobh — c* relation, 
and our relatively short binary decay timescales. By con- 
trast, our predicted number counts are about an order of 
magnitude smaller than given by Wyithe & Loeb (2003). 
This is because they do not consider the shift in frequency 
during the inspiraling phase, and compute the counts at a 
fixed frequency of 10 -3 Hz. Their simplified treatment has 
the effect of artificially boosting (by as much as a factor 
of ~ 50) the number of detectable events at z > 7. 

Lastly, we were concerned that, in the LISA window, the 
GWB we computed is entirely due to GWs emitted at the 
last stable orbit of the coalescing pair, a regime the New- 
tonian approximation may fail to describe adequately. To 
quantify this we have implemented a post-Newtonian (or- 
der 2) approximation describing the GW emission from bi- 
nary systems (Blanchet 2001). We find that differences in 
the spectrum of a coalescing MBH binary are quite small, 
of order 10%, and the overall results are, basically, unaf- 
fected. 

5. SUMMARY AND CONCLUSIONS 

We have computed the gravitational wave signal (in 
terms of the characteristic strain spectrum) from the cos- 



mological population of inspiraling MBH binaries pre- 
dicted to form at the center of galaxies in a hierarchical 
structure formation scenario. The assembly and growth of 
MBHs was followed using Monte Carlo realizations of the 
halo merger hierarchy from the present epoch to very high 
rcdshifts, coupled with semi-analytical recipes treating the 
dynamics of the (inevitably forming) MBH binaries. We 
find that the broad band GWB spectrum can be divided 
into three different regimes. For frequencies < 10~ 10 Hz, 
the GWB is shaped by MBH binaries in the gravitational 
slingshot regime, i.e. the orbital decay is driven by scat- 
tering off background stars. In the intermediate band, 
10~ 9 < / < 10~ 6 Hz, GW emission itself dominates poten- 
tial energy losses and the strain has the "standard" /~ 2 / 3 
behaviour. Finally, for / > 10~ 6 , the background signal is 
formed by the convolution of the emission at the last stable 
orbit from individual binaries. In the first two regimes, the 
strain is dominated by low redshift events (z < 2) involv- 
ing genuinely supermassive holes. At larger frequencies 
the main contribution comes from lighter and lighter bi- 
naries. Such background should dominate that predicted 
from a population of cxtragalactic white dwarf binaries. 
As already pointed out by Wyithe & Loeb (2003), in the 
nHz regime probed by pulsar timing experiments, the am- 
plitude of the characteristic strain from coalescing MBH 
binaries is close to current experimental limits (Lommen 
2002). 

In the LISA window (10~ 6 < / < 0.1 Hz), the main 
sources of GWs are MBH binaries in the mass range 
10 3 <mi < 10 7 M . With a plausible lifetime of ~ 3 
years, LISA will resolve the GWB into s=s 60 discrete 
sources above S/N = 5 confidence level. Among these, 
w 35 are bursts, i.e., binaries caught in their final inspiral- 
ing phase. We predict that most of the observable events 
will be at redshifts 2 < z < 7, with only a handful at larger 
rcdshifts. While LISA will make it possible to probe the 
coalescence of early black hole binaries in the universe, it 
may not be able to observe the formation epochs of the 
first MBHs. 

We have benefitted from several discussions with V. 
Gorini. Support for this work was provided by NASA 
grant NAG5-11513 and NSF grant AST-0205738 (P.M.). 
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